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Recently, Yin et al. [Eur. Phys. J. B 49, 205 (2006)] introduced an efficient small-world net- 
work traffic model using preferential next-nearest neighbor routing strategy with the so-called path 
iteration avoidance (PIA) rule to study the jamming transition of internet. Here we study their 
model without PIA rule by a mean-field analysis which carefully divides the message packets into 
two types. Then, we argue that our mean-field analysis is also applicable in the presence of PIA 
rule in the limit of a large number of nodes in the network. Our analysis gives an explicit expression 
of the critical packet injection rate Rc as a function of a bias parameter of the routing strategy a in 
their model with or without PIA rule. In particular, we predict a sudden change in Rc at a certain 
value of a. These predictions agree quite well with our extensive computer simulations. 
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I. INTRODUCTION 



Complex networks with small-world property exist in 
many natural and social systems, such as food web, the 
internet the world wide web 3], and the world-wide 

airport network (WAN) [4]. In 1999, Barabasi and Albert 
proposed a scale-free growing model (BA network) with 
a preferential attachment mechanism to mimic a growing 
small- world network in the real world Their model 
stimulated the interest of the physics community to study 
complex networks by statistical physical means @, 
One of the goals of these studies is to understand the 
dynamical processes taking place behind the underlying 
structure. 

It is instructive to study the traffic capacity of a net- 
work. We may start by considering the simple-minded 
situation in which message packets are injected randomly 
into the nodes of the network at a fixed rate. Each 
packet has a randomly assigned destination node. And 
each node in the network has a finite message-forwarding 
rate. Clearly, an important factor affecting network traf- 
fic capacity is the routing strategy, namely, how each 
node forwards its out-going message packets to its nearest 
neighbors. The performance indicator is the maximum 
free-flowing traffic capacity characterized by the critical 
packet generation rate Rc- More precisely, Rc is the 
supremum number of new packets that can be injected 
into the network per unit time step without causing con- 
gestion [1, Q. Here, congestion means that the average 
rate of change of the number of packets in some node 
is positive. (Actually, this performance indicator is not 
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overly stringent for the model investigated in this paper 
as we find that the number of packet in almost all nodes 
steadily increases over time without saturation whenever 
R > Rc-) The more efficient the routing strategy, the 
larger the value of Rc- For a sufficiently large random 
network, the routing strategy cannot depend on the net- 
work topology because this information is not available 
to each node. Thus, it is reasonable to confine ourselves 
to study local routing strategies. 

Perhaps the simplest local routing strategies are the 
ones that use information on the nearest neighbors of in- 
dividual node [ll| . Recently, based on these nearest- 
neighbors-based strategies, a new routing strategy called 
the preferential next-nearest-neighbor (PNNN) searching 
strategy was proposed by Yin et al. [l^ in which the 
performance is better than those using nearest neighbor 
routing. As the name suggests, in PNNN, a message 
packet looks for its destination among the next nearest 
neighbors of the node it currently stays. If the desti- 
nation cannot be found in this way, the message packet 
will be forwarded to a neighboring node by a biased ran- 
dom walk with a preferential probability which depends 
on a parameter called preferential delivering exponent a. 
To speed up packet delivery, Yin et al. added in their 
routing strategy the path iteration avoidance (PIA) rule, 
which states that a packet cannot travel through an edge 
more than twice. 

As a model of scale-free network traffic with potential 
applications in the internet and the world wide web, the 
use of PIA rule is problematic. A message packet, un- 
like a human driver, cannot automatically remember the 
path it has traveled. This additional piece of information, 
whose length grows linearly with the time since creation 
of the packet, may either be stored in, say, a central reg- 
istry, or attached to the message packet itself. Thus, 
the cost of inquiring this information from the registry 
or transmitting it through an edge alongside with the 
message packet cannot be ignored. Furthermore, addi- 
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tional computational cost, which also scales linearly with 
the time since the creation of the packet, is needed for 
a node to process this historical path information in ac- 
cordance with the PIA rule. All these factors make the 
effective message packet forwarding rate a function of 
the time since the packet creation. Unfortunately, the 
PNNN routing strategy of Yin et al. f\3\ does not take 
these extra communication and computational costs into 
account. This is why we believe that PIA rule is not very 
realistic. 

In Sec. ini we briefly review the network traffic model 
proposed by Yin et al. [l^] using the PNNN strategy. 
Then we perform mean-field analytical calculation for the 
dynamics of their model with and without the PIA rule 
in Sec, mil In both cases, we find an abrupt change in the 
dependence of Rc on a at certain value of a. We also give 
the physical reason behind such change. In Sec. IIVI we 
compare the mean-field calculations with our extensive 
numerical simulation results of Rc against a. We also 
show in this Section that the network size used in Yin et 
aL's numerical simulations is not large enough to reveal 
the thermodynamic behavior of their model. Finally, we 
give a brief summary and discuss the effectiveness of the 
PNNN strategy in Sec. El 



II. THE PNNN+PIA AND PNNN-PIA MODELS 

Yin et al. proposed and studied the following network 
traffic model on a BA network J^]. (Here we call their 
model with and without the PIA rule PNNN+PIA and 
PNNN-PIA, respectively.) Their model consists of a ran- 
dom but fixed BA network with N nodes. We denote the 
set of all nodes in this network by V. We further denote 
the degree of the node having the least (greatest) number 
of nearest neighbors in the network by k^i^x (fcmax)- That 
is to say. 



and 



max ki 



(1) 



(2) 



where ki is the degree of the node i. Recall that BA 
network is generated by connecting each newly added 
node to m existing nodes in a careful way Hence, 



k ■ — 



(3) 



Further recall that during the generation of a BA net- 
work, the average degree of the node added to the net- 
work toiapsc ago equals m(iV/iciapsc)"^/^ W\- Thus, 



(4) 
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We denote the adjacency matrix of the network by A. 



In PNNN-I-PIA and PNNN-PIA, each node has an un- 
limited buffer, known as load, to store packets. At each 
time step, each of the R packets is added to a randomly 
chosen source node of the network with a randomly cho- 
sen destination node. Note that simulations reported by 
Yin et al. in Ref . [12] were performed by considering only 
integer values of R. In contrast, we allow a real-valued 
R. More precisely, we inject a message packet into a node 
with probability R/N in each time step. 

Each node can send out at most C > 1 packets to its 
nearest neighbors using the first-in-first-out rule. That 
is to say, packets entering a node first will be sent out 
first. Each out-going packet first searches through all the 
next nearest neighbors of the node to which it currently 
belongs. If its destination is located in this search, the 
packet will be forwarded to one of the neighbors connect- 
ing the destination and the current node. And in the next 
time step, this packet will be forwarded to the destination 
and then removed from the network. If the destination 
of an out-going message packet cannot be found in such 
a search, it will be randomly forwarded from its current 
node (say, node i) to one of the neighbors (say, node j) 
with probability 



H„- = 



kf 



(5) 



That is, Aij 
nodes i and j. 



1(0) if there is an (no) edge between 



where a is a fixed parameter known as the preferential 
delivering exponent. Note that the sum in the above 
equation can be regarded as a restricted sum over the 
nearest neighbors of the packet's current node i. 

The only difference between PNNN-HPIA and PNNN- 
PIA is that PIA rule is present in the former model while 
absent in the latter. Recall that PIA rule demands each 
packet to travel through the same edge at most twice [l^ ■ 
In the event that a message packet has nowhere to go 
due to the PIA rule, the packet will be removed from 
the network. And for a G [—4,2], only a very small 
percenta ge o f packets are removed from the network in 
this way [13| . 

Clearly, historical path information of a packet is 
needed to decide where it will go in the next time step 
with the adoption of PIA rule. As we have mentioned 
in Sec. [H extra communication and processing costs are 
required to forward a message packet together with its 
historical path information in the network to its neigh- 
boring node. Thus, it is less efficient to forward an old 
packet than a newly created one. In this respect, PIA 
rule is not consistent with the rule that the message for- 
warding capability of a node is independent of the age of 
the forwarding packets. This is a serious problem because 
Yin et al. found by numerical simulation that the packet 
lifetime, which is the time between its injection and re- 
moval, roughly obeys a power law distribution 112)]. 

Although Yin et al. has briefly studied the PNNN- 
PIA model numerically in Ref. [l^, their focus was on 
the PNNN-fPIA model. They found that the critical 
packet generation rate Rc is increased by adopting the 



3 



PIA rule. More importantly, using numerical simula- 
tion up to TV « 5000 with R restricted to integers only, 
they found that Rc is a decreasing function of a for the 
PNNN+PIA model. In addition, based on their sim- 
ulations in the range a S [—4,2], they believed that 
for a fixed N, the value of Rc is a constant whenever 
a < -2 [H. 

An interesting common feature of the PNNN-I-PIA and 
PNNN-PIA models is that as long as there are no more 
than C message packets staying in a node at any time, 
the message packets behave like independent particles 
in the sense that their motions in the BA network are 
independent of each other. This property is important in 
our subsequent discussions. 



III. MEAN-FIELD ANALYSIS 
A. The PNNN-PIA Model 

We try to calculate the Rc against a curve for the 
PNNN-PIA model by mean-field approximation. The 
validity of the approximations made in our calculation 
will be discussed and justified in Sec. IIVI Let us begin 
by classifying the packets into two types. A packet is 
called a destination located packet (DLP) if it has suc- 
cessfully found a path to its destination. By the rules of 
PNNN-PIA, the destination of a DLP must be one of the 
nearest or next nearest neighboring nodes of its current 
location. Otherwise, the packet is known as a destina- 
tion seeking packet (DSP). Since a newly injected packet 
has not found out its path to the destination yet, it must 
a DSP. A DSP moves randomly to its neighboring node 
with probability given by Eq. We denote the num- 
bers of DLPs and DSPs in node i at time t by ni^i{t) and 
ns,i{t), respectively. 

We say that a network is in free-flow state if each node, 
on average, can forward all its loads in the next time step. 
(In other words, the average load of each node is at most 
C in each time step.) In this case, node i can, on average, 
send out all its ns^i{t) message packets at any time t. At 
the same time, node i receives, on average, R/N DSPs by 
packet generation. Since the number of 4-cycles in a BA 
network scales like [m\og{N)/2]'^ /4 yu,], the probability 
that two next nearest neighboring nodes are connected 
to more than one common node goes to in the large N 
limit. So the number of next nearest neighbors for node i 
is approximately equal to J2jev ^iji^j ~ 1) large 
N limit. Consider a DSP that reaches the node i for the 
first time. Then, the probability $i that it can locate a 
path to its destination in the next time step is given by 



X! 

j<£V 



(6) 



In contrast, suppose the DSP has reached the node i 
more than once, then it has no chance to find the path 
to its destination in the next time step as the next nearest 
neighbors of node i has been searched during its previous 
visit to node i. Let Xj be the average number of visit of 
a DSP to node j given that it has visited node j at least 
once. Then, using the mean-field approximation similar 
to that used in Ref. [15,], the number of DSP in free-flow 
state satisfies 



Jt 



R 

N 



(7) 

Note that 1 — ^j/Xj is the probability that a DSP at 
node j will not change to a DSP in the next time step. 
Thus, the last term in the R.H.S. of the above equation 
is the average number of DSPs received by node i from 
its neighbors. 

We want to study the equilibrated distribution of DSPs 
for a typical node in the free-flow state as a function of 
the degree of the node. And we do so by investigating 



ns{k) 



(8) 



where 6ki,k is the Kronecker delta and (•■■)( represents 
the time average of its argument. Upon equilibration. 



R 

N 



^A,, (n,j(i))^n,, 1 



(9) 



Although BA network does not show assortative mix- 
ing 16], it exhibits non-trivial but weak degree-degree 
correlation between neighboring nodes @. Combined 
with the fact that the trajectory of a DSP is history in- 
dependent for PNNN-PIA, it makes sense to ignore this 
degree-degree correlation in our mean-field analysis. By 
ignoring this correlation, we know that for any function 
f{k), 



(10) 



where I? is a functional of /. Most importantly, V is 
independent of kj. As a result, using Eqs. (O-®, we 
can re-express Eq. ^ as 
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risik) 



R 

N 



E 



Ski,kAijn-s{kj) 



Efev^J^(^^ - 1) 



R 

N 



Dk 



a+l 



(11) 



r 



for some D > Q independent of k. Of course, D depends 
on a and N . 

To derive the mean field equation for ni^i{t) in the free- 
flow state, we consider a DSP currently located at j, 
which is a neighboring node of i. Suppose this is the 
first time for this DSP to visit a neighboring node of i. 
Then, on average, the chance for this packet to turn into 
a DLP and then forwarded to node i in the next time 
step equals ^j[{ki - l)/J2i,Aje{kt - 1)] = (h ~ l)/N 
in the large N limit. In contrast, if this is not the first 
time for the DSP to visit a neighboring node of i, then 
it has no chance to be forwarded to node i as a DLP in 
the next time step. This is because the DSP should have 
converted into a DLP after its first visit to a neighboring 
node of i. Suppose a DSP was located at a neighboring 
node of i at time step < — 1. Suppose further that this 
packet is forwarded to node i at time step t. Then it 
must be found in a neighboring node of i at time step 
t + 1. In contrast, suppose the packet is forwarded to a 
node other than i at time step t, then the chance that it 
will come back to a neighboring site of i will be roughly 
proportional to fc^. Hence, the average number of times 
for a DSP to visit neighboring nodes of i given that it 
has visited a neighboring node of i once equals (/i + I'ki) 
for some /i > 1 and > independent of ki. 

It is obvious that fj, is independent oi N. In what 
follows, we argue that v is also independent of A^. BA 
network is a small world network without showing any 
assortative mixing 16]. And the packet forwarding rule 
in Eq. ([5]) does not depend on the historical path of the 
packet. So, message packets are essentially performing 
random walk in the network in the first few steps after 
its injection. Consequently, the probability distribution 
of the first return time of a random walker should scale 



like for some ^ > and sufficiently small t [17[. More- 
over, ^ is independent of N. On the other hand, when 
t is approximately greater than the average square dis- 
tance between two nodes in the network {(P), finite size 
effect of the network will affect the probability distribu- 
tion of the first return time of a random walker so that 
the scaling will no longer be valid. Indeed, this is 
what Almaas et al. have found in their numerical study 
of the first return time for random walk in a certain small 
world network. More importantly, they found that the 
probability distribution of the first return time collapses 
to a single scaling relation by rescaling both the first re- 
turn time t and the probability P by ((P) [11]. Since the 
packet lifetime r scales roughly as (cP), we conclude that 
V is independent of N. Nevertheless, both ^ and v are 
functions of a. But the form of Eq. (O assures that /i 
and ly are not sensitively dependent on a in the sense that 



fi and v scale polynomially instead of, say, exponentially 
with a. 

Utilizing all these information, we may write the mean 
field equation for ni^i{t) in free-flow state as follow: 



dt 



N{fi + vki) 



Note that the average number of DLPs for a typical 
degree k node in the free-flow state upon equilibration: 



ni{k) 

satisfies 

ni{k)^5u^^k 



SiGV ^ki,k 



(13) 



fc- 1 



+ vk) 



Sk.MAjrisikj). (14) 



By ignoring the degree-degree correlation between neigh- 
boring nodes as in the derivation of the scaling relation 
for ns{k), we have 



ni{k) 



(ns(fc»)),gvfc(fc- 1) 
N{p + uk) 

("-s(fci))igv k 



(15) 



(16) 



As we shall see in Sec. IIVI the value of v is of order of 
0.01 for most values of a. Thus, for network size N < 
5000 such as those used in the simulations reported in 
Ref. [T^, ni{k) varies quadratically rather than linearly 
in most of the domain Jfcmin , ^max] • In this respect, Yin 
et aL's numerical results did not reflect the properties of 
the system in the large N limit. We shall discuss more 
along this line in Sec. IIVI 

Upon equilibration, the average number of packet re- 
siding on a typical degree k node equals 



n{k) = ns{k)+ni{k) 



R 

N 



Dk 



a + l 



{ns{ki))^^^k{k- 1) 



(17) 



N{/j + uk) 

^e^- (18) 



Nv 



in the large N limit. 



1. A Simplifying Assumption 

In this Subsection, we make the simplifying as- 
sumption that the expressions of ns{k) and ni{k) in 



5 



Eqs. (fTTj) and ((T5)) are exact throughout the entire do- 
main [fcininj ^max]- Then, it is clear that Eq. pT]) is an 
increasing function of k for a > — 1. Hence, the max- 
imum value for the last line of Eq. psp in this domain 
is attained when k = fcmax- And in the case of a < — 1, 
Eq. (jl7p is a continuous function with one local minimum 
point in the interval [fcmin, ^max]- So, again in this inter- 
val, n{k) attains its maximum value at the boundary. To 
find out the exact location at which the maximum value 
is attained, we have to find an expression for {ns{ki)) 
first. 

According to Albert and Barabasi, the probability dis- 
tribution of nodes of degree k for a BA network is given 



by 



(19) 



where E is the normalization constant and 7 = 3 [6| . So 
the normalization constant E can be rewritten as 



(7-l)m'^-i. (20) 



Using our assumption that Eq. (jlip is valid over the en- 
tire interval [fcmin , ^max] , we arrive at 




(n,(fc,)).ev= / pik)nMdk^ a-'y + 2 ^^^^ 



fcmin 



in the large N limit provided that a 7^ 7 — 2 = 1. 

By substituting Eqs. ([3]), 11]) and pTjl into Eq. (fT8|) together with the fact that v is independent of N and is not 
sensitively dependent on a, we find 



4-1 f r 2m ri - iV("-i)/21 (1 - iVi/2) 

7Vz/(l — a) 



> (22) 



Thus, the max- 
n provided that 



in the large N limit whenever a < — 1. 
imum of n(fc) is attained at fc = k^ 
a < —1 and TV — > 00. 

To summarize, the maximum value of n{k) is always 
attained either at A: = fcmin or fc = fcmax- By denoting 
the value fc at which n(fc) reaches its maximum value by 
fcc, we have 



lim fcc 



fcmax if a > -1, 

fcmin otherwise. 



(23) 



And from Eq. ([22)l . for any fixed iV > 0, there is a critical 
value of a = ttc < — 1 above (below) which fcc = fcmax 



(fcc 



Besides, 



lim 

N^oc 



= -1. 



(24) 



There is an important consequence of the above find- 
ings. By gradually increasing the packet injection rate i?, 
the first congested node must be the one with the largest 
value of n(fc). Therefore, the critical packet injection 
rate Rc is reached when congestion occurs at a smallest 
(largest) degree node whenever a < ac (a > etc). This 
change in the type of node that congests first upon a 
gradual increase in R results in the discontinuity of Rc 



at a = ckc- Actually, Rc attains its maximum value at 
a — ac in the large N limit. To see why, we consider the 
situation when the network is at its maximal capacity. 
In this situation, R — Rc and the maximum number of 
packets in some nodes should be C. From Eq. ^TT} . Rc 
satisfies 



C 



D{l l)m^-i (k^-J+^ k^-'^') fcc(fcc - 1) 



N{a~j + 2){fi + iykc) 

a + l 



+ 



(25) 



We need to eliminate D in order to simplify the above 
equation. We proceed by considering the average number 
of packets reaching their destinations in each time step at 
equilibrium. This number is equal to the average number 
of packets injected into the system at each time step. 
Therefore, 



R 



Np{k)ni{k)dk. 



(26) 



From Eqs. ©-(U), US]) and (HHl-lIlT]), we know that the 
critical packet injection rate equals 



6 



Rc 



E{ns{h)).^^ 



k - 1 



/i + 



In 



,7n J 



a - 7 + 2 

in the large A'' limit provided that a ^ 7 — 2 = 1 . 

By using Eq. ([27| to eliminate Z? in Eq. ((25|l . we find that a sufficiently large N and a 7^ 7 — 2 = 1, 



Rn 



CN{j - l)2m"+T (^N 2* - ij 




(7 - l)2m°+T 


- 1) 


S + (a - 7 + 2)iVfc?+i + (7 - 


(iV^ - 1) 


fce(fec-l) 

fi-\-i'kc 



mm 

fce{m,-mVlV} I 



(27) 



(7 - l)2m"+')' (iV^^^ - 1) S + (a - 7 + 2)iVfc"+i + (7 - l)m"+i (. 



k{k-l) 



m 



where 



than simply substituting a = 7 — 2 into Eq. 



■In 



(/i + vm)\/N 



N 



N 



— ■ (29) 
/im 



The above equation is not only valid for the generic case. 
It is straight-forward to go through the same derivation 
to show that Eq. (pS)) is also valid for the singular case of 
a = 7 — 2 as long as we take the limit a — )■ 7 — 2 rather 



Although the functional forms of /i and i' are not easy 
to determine, the facts that they are independent of 
and are not sensitively dependent on a are already suffi- 
cient for us to make the following remark on the general 
trend of Rc- For sufficiently large N, R^ is an increas- 
ing (decreasing) function whenever a < ac {a > ac). 
Besides, 



lim Rc = 



7.7-1 



■In I 



1 



> for a < -1, 



(30a) 



lim Rc = lim 

iV— >oo N^oo 



C(7- 1)2to 



^ln( 



1 



(7 



2)7V("+i)/2 



for - 1 < a < 1, 



(30b) 



and 



lim Rc = lim 



C(7-l)2m7-i ii+fl In (ii±i^) 



1 



{a-j + 2)N 



= for a > 1. 



(30c) 



In other words, in the thermodynamic limit, the change 
in the type of nodes that is congested first results in the 
maximum point of the a — Rc curve at a = ac- 

We may understand the occurrence of this maximum 
point as follows. Clearly, there are much more small 
degree nodes than large degree ones in a BA network. 
Since all nodes of different degree have the same message- 
forwarding capability, one may attempt to increase Rc by 
preferentially forwarding the DSPs to small degree nodes 



by setting a<0. If ac < a < 0, the bias towards sending 
DSPs to small degree nodes is not yet sufficient. Hence, 
jamming at R = Rc occurs in the largest degree node 
because too many DLPs move to this node per unit time 
step. In contrast, if a < ac, the bias towards sending 
DSPs to small degree nodes is too strong that the small- 
est degree nodes are jammed by the infiux of DSPs. In 
this respect, it is not surprising for our mean-field cal- 
culations to find that Rc is an increasing (decreasing) 
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FIG. 1: Log-log plot of the distribution of number of DSPs 
ris against the degree of node k in black dots for PNNN-PIA 
with m = 4, C = 1 and R = 5 for network size A'' = 10^ and 
preferential delivering exponents a. The dashed line with 
slope a + 1 in each subplot is drawn for comparison purpose. 



function of a < (a > etc). 

They may break down near A^min and /^niax- a re- 
sult, the expression for D in Eq. (|27|) should only be 
regarded as a trend indicator. Besides, upon gradual in- 
crease in the packet generation rate R, the first node 
to be congested may no longer be the one whose degree 
is fcmin or fc,nax- Nevertheless, the maximum point on 
the a — Rc curve due to the change of the kind of node 
that is congested first is robust and generic as it is stable 
upon small change in n{k). Of course, the expression for 
Rc in Eq. (|28p and the value of ac will be affected as 
a consequence of the break down of Eqs. pT|) and 
Fortunately, as fj, and v are independent of N and are 
not sensitively dependent on a, we conclude that S is 
almost N independent in the large N limit. More im- 
portantly, within about 10% accuracy, we may regard 
S as independent of a. Thus, the general trend of Rc 
expressed in Eqs. (|30ap - p0cp is still valid. That is to 
say, for sufficiently large N, Rc is approximately propor- 
tional to 1/(7 — a — 2) = 1/(1 — a) for a < ac- And 
the proportionality constant is independent of N. Be- 
sides, Rc ~ 1/[(1 - a)7V("+i)/2] for -1 < a < 1 and 
Rc ~ l/[N{a — 1)] for a > 1. We are going to test 
these predictions using large scale numerical simulations 
in Sec. |TVl 



B. Implications To The PNNN+PIA Model 

1. Beyond The Simplifying Assumption 

In reality, Eqs. pT|) and are not exact. Although 
it is much harder to modify the mean-field analysis in 
Sec, nil Al to take the PIA rule into account, we can still 



FIG. 2: Log-log plot of the distribution of number of DSPs Us 
against the degree of node k for PNNN+PIA. All parameters 
used are the same as those in Fig. [T] 



argued the behavior of the PNNN+PIA model qualita- 
tively. First, we consider the effect of PIA rule on 7is(fc). 
Clearly, PIA rule makes Hji in Eq. ([7|) historical path 
dependent. Thus, we can no longer apply the trick in 
Eq. pO)) to give a simple expression for ns(k). Neverthe- 
less, we may argue the behavior of ns{k) as follows. In 
the case of a < 0, Eq. ^ implies that packets are pref- 
erentially being forwarded to small degree nodes. How- 
ever, the PIA rule forbids a packet to travel through the 
same edge more than twice. Therefore, compared with 
the situation without the PIA rule, a packet is less likely 
to be forwarded to a small degree node on average. On 
the other hand, the PIA rule has relatively little effect 
on high degree nodes. This is because of two reasons: 
first, packets are less likely to travel to these nodes; and 
second, packets located at these nodes generally have a 
large number of possible nodes to be forwarded to in the 
next time step. Thus, we expect that the same scaling 
behavior for ns{k) found in Eq. pTjl is observed in the 
presence of PIA rule. Nonetheless, the domain of k in 
which this scaling law holds is reduced as the value of 
the lower cutoff of the scaling law increases as a conse- 
quence of the PIA rule. Furthermore, below this lower 
cutoff point, the value of ns{k) is smaller that the case 
when the PIA rule is not adopted. 

Applying similar arguments in the previous paragraph 
to the case of a > 0, we conclude that it is more likely to 
forward a packet between two large degree nodes. Since 
the number of nodes with degree k decreases as k in- 
creases, the combination of Eq. ^ and the PIA rule 
will decrease (increase) the value of ns(fc) in Eq. (jlip for 
k < fcniax [k ^ Aimax and fc ^ fcmin)- Therefore, the do- 
main in which the scaling behavior of Eq. pd|) holds only 
for k <C fcmax and k ^ fcmin- To summarize, we have ar- 
gued the validity of Eq. ([TT|) in the large TV limit for the 
PNNN+PIA model over a reduced domain of fc. In ad- 
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FIG. 3: Log-log plot of the distribution of number of DLPs 
ni against degree of nodes k for PNNN-PIA. The solid curve 
in each subplot is the prediction according to Eq. (I15|l with ^ 
and u treated as free fitting parameters. And the dotted line 
in each subplot is the asymptote of the corresponding solid 
curve. All parameters used in the simulations are the same 
as those in Fig. [1] 



dition, the value of {ns(ki))i^Y decreases in the presence 
of PI A rule. 

How about the effect of PIA rule on ni{k)7 The PIA 
rule surely reduces both /i and by forbidding exces- 
sive routing through the same edge. Besides, the value 
of (ns(fci))igv is also reduced. But interestingly, unlike 
Eq. ([7]), the presence of PIA rule in no way affects the 
functional form of Eq. pS]) as the derivation of Eq. ^2]) 
is also valid in this case. This is because the PIA rule 
cannot prevent a DLP from reaching its destination un- 
less the distance between the destination and the initial 
generation point of the packet is less than two. (This is 
because at the first instance when a DSP is forwarded 
to a node i with distance two from the packet destina- 
tion. The DSP will turn into a DLP in the next time 
step. More importantly, this message packet must never 
pass through any shortest path connecting node i and the 
packet destination. Hence, the PIA rule does not prevent 
this packet from moving along this shortest path.) And 



FIG. 4: Log-log plot of the distribution of number of DLPs 
ni against degree of nodes k for PNNN-I-PIA. The detailed 
procedure is adapted from the descriptions in Fig. |3l 

the probability for such case is negligible in the large 
N limit. Note that ni{k) is more seriously affected by 
{ns{ki))ieY than by /i or v. So, we expect that ni{k) 
decreases with the introduction of PIA rule. But its per- 
centage decrease is not as large as that of ns{k^ax)- 

We now move on to study the effect of PIA rule on 
the values of ac and Rc- Recall that without PIA rule, 
ttc — —1. Let us consider the case of a > — 1 first. In this 
case, both ns{k) and ni{k) are increasing functions of k 
for a > — 1 with or without PIA rule. So, upon a gradual 
increase in the packet injection rate, the first node to 
congest must be the one with a large degree. From the 
arguments in this Subsection, we know that for k w fcmax, 
n(k) = ns{k) + ni(k) decreases with the introduction of 
PIA rule. Hence, the critical packet injection rate Rc 
increases with the introduction of PIA rule. Certainly, 
the percentage increase in Rc depends on the values of 
N , m and a used; and the above arguments in no way 
imply that the percentage change is huge. Indeed, it is 
quite possible that the increase in Rc is negligible in some 
cases. 

In contrast, when a < —1, Eq. (fTB|) together with 
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FIG. 5: Plots of fi against a for (a) PNNN-PIA and 
(b) PNNN+PIA for various values of N. 
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(b) PNNN-I-PIA for various values of N. 



IV. COMPARISON WITH OUR NUMERICAL 
SIMULATIONS 



the arguments in this Subsection tell us that n(fcniin) ~ 
nsikmin) decreases more rapidly than n(fcmax) ~ 
«K^max) ~ ('^s(fci))iev ''^max with the introduction of PIA 
rule. Consequently, Rc increases in the presence of PIA 
rule. More importantly, for a finite N, one may find an 
a slightly less than —1 such that n(fcniin) < '^(fcmax)- In 
other words, a large degree instead of a small degree node 
gets congested at Rc for this value of a. Therefore, we 
conclude that ac decreases and Rc increases in the pres- 
ence of PIA rule. Note that once again the decrease in 
ttc may be insignificant in some cases. 

Finally, we expect that the general trend of Rc for 
PNNN-PIA described in Sec. HUB II also applies to 
PNNN-I-PIA. Obviously, our predictions are different 
from the numerical results of Yin et al. reported in 
Ref. which claimed that Rc was a decreasing func- 
tion of a for PNNN-hPIA. In Sec. Ml we show that this 
is partly due to the fact that the network size N used 
in their simulation is not large enough so that finite-size 
effect seriously affects their conclusions. 



We want to check the validity of our mean-field analysis 
reported in the previous Section as well as to to under- 
stand the origin of the discrepancy between our present 
work and the numerical results obtained by Yin et al. in 
Ref. [l^ . And we do so by performing numerical simula- 
tions using larger values of N. Moreover, unlike Ref. 
we allow R to take on non-integer values. 

Perhaps one of the reasons why Yin et al. reported 
numerical simulations of PNNN-I-PIA up to iV = 5000 
only is that a lot of memory is needed to store the 
message packets present in the network as well as their 
historical paths. In fact, this straight-forward numerical 
simulation method is not practical for N > 10000. Here 
we introduce a much less memory intensive way to nu- 
merically find Rc. Observe that the connectedness of BA 
network and the message forwarding rules of PNNNiPIA 
make the message packets in PNNNiPIA crgodic. Also, 
recall from Sec. |TT] that message packets behave like in- 
dependent particles as long as there are no more than 
C packets staying in a node at any time. Although oc- 
casionally more than C message packets may be present 
in a node in the free-flow phase, by ergodicity we expect 
that the statistical properties of PNNNiPIA below the 
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FIG. 7: The average number of packets n against degree of 
nodes k for PNNN-PIA at R — Re- All parameters used in 
the simulations are the same as those in Fig. [1] 



FIG. 8: The average number of packets n against degree of 
nodes k for PNNN+PIA at R ^ Rc. All parameters used in 
the simulations are the same as those in Fig. [1] 



critical packet injection rate Rc can still be simulated 
by regarding each message packet as independent par- 
ticle throughout. Therefore, the statistical behavior of 
PNNNiPIA for R < Rc can be found as follows: We first 
numerically simulate the ensemble-averaged time evolu- 
tion of a particular free-flow phase situation in which 
there is exactly one message packet in the network at all 
times. By ergodicity, the ensemble-averaged number of 
packet present in a node obtained in the above simula- 
tion equals the (time-averaged) number of packet in that 
node when the packet injection rate i? is 1/(t), where (r) 
denotes the mean packet lifetime. (This choice of R does 
not contradict with the prediction of Eqs. (j30bp and (|30c[) 
that i?c in the limit of large N whenever a > ac- 
This is because the mean packet lifetime (t) scales like 
with P > 2.) Below the critical packet injection 
rate Rc, the distributions ns{k), ni{k) and n{k) are di- 
rectly proportional to the packet injection rate R. Conse- 
quently, Rc is equal to C/{t) maxfc n{k) where max/j n{k) 
is the maximum value of n(k) over all k for the case of 
R = 1/(t). Clearly, this method can compute Rc accu- 
rately and efficiently. As only one message packet is used 
at any time in the simulation, this method requires much 
less memory than the straight-forward numerical simu- 



lation approach. We further verify the validity of this 
ensemble-averaged simulation method by successfully re- 
producing the numerical simulation results reported by 
Yin et al. in Ref. [12] (modulo the fact that they re- 
stricted R to integers). (Actually, the value of m used 
for their PNNN-PIA simulation is 4 instead of 5 19'].) 
Therefore, we adopt this new method in our subsequent 
numerical studies. 

While the simulations of Yin et al. in Ref. [13 was 
performed in for a e [—4,2], ours is done in a slightly 
large parameter range of [—4,4]. Actually, we find that 
the bias in forwarding a DSP according to Eq. ([5]) for ]a] 
close to 4 is already so high that the data obtained from 
our simulations are no longer very reliable. And reliable 
results for \a\ > 4 has to be obtained by much longer 
simulation time with the aid of a higher precision pseudo 
random number generator. 

Let us begin by checking the validity of our assump- 
tions made in Sec. IIIII Figs. [T] and [2] show typical 
ns{k) curves obtained from our numerical simulations of 
PNNN-PIA and PNNN+PIA, respectively. They show 
that risik) indeed follows a power law with exponent a+l 
over most of the parameter range for PNNNiPIA for suf- 
ficiently large N . Furthermore, the domain of validity of 
the power law is reduced with the introduction of PIA 
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a a 



FIG. 9: The Rc against a curve for PNNN-PIA with C = 1 
and (a) m = 4, and (b) m — 5. The dashed curve in each sub- 
plot is our mean field analytical prediction based on Eq. (|30a|) . 
More precisely, the dashed curve is the best fit curve obtained 
fi-om Eq. pOa|l by treating H as a free parameter independent 
of Q. 



rule. More importantly, in the case of PNNN+PIA, the 
ways how ns{k) deviates from the power law for small and 
large k are consistent with our predictions in Sec. IIIIBI 
That is to say, ns{k) is less (greater) than the value ob- 
tained by Eq. pT|) for k « /cmin {k ~ fcmax)- In this 
respect, our assumption of ignoring degree-degree corre- 
lation between neighboring nodes in obtaining ns{k) is 
not bad. 

Next, we examine the validity of Eq. ^T5\i for 
PNNNzbPIA. Figs. [3] and H plot ni as a function of 
k obtained from our simulation of PNNN-PIA and 
PNNN+PIA, respectively. Our simulation results for 
ni{k) agree quite well with the solid curves, namely, our 
mean field prediction given by Eq. (fT5| . The dotted lines 
in Figs. [3] and d] show the asymptotic behavior of the 
solid curve in the limit of large k. By comparing our 
simulated data points with the dotted lines, we find that 
for N as small as 1000, ni{k) does not reach the linear 
scaling regime at all. And for N = 100000, n;(fc) attains 
linear scaling for k > 200. In fact, we discover from our 
simulation that ni{k) scales like a linear function of k 



FIG. 10: The Rc against a curve for PNNN+PIA. Parameters 
used are the same as those in Fig. [Q] 



around k < fcmax only when N > 10000. 

As shown in Figs. [5] and [51 /i and v are independent of 
N for PNNN+PIA provided that TV > 7000 and |a| < 3. 
We believe that the discrepancy for /z when N — 1000 in 
Fig. [5] is the result of finite size effect. And as we have 
already discussed earlier in this Section, we think that 
the discrepancies for fi and v for la] > 3 are due to the 
limitations of our simulation time and pseudo random 
number generator used. In any case. Figs. [5] and [S] verify 
that n and v are not sensitively dependent on a. In fact, 
/i and v are of order of 1 and 0.01 respectively over most 
of the range of a we have studied. And in line with our 
expectation, /i and v decrease with the introduction of 
PIA rule. 

Figs. [7] and [5] depict the general trend of ris(fc), n;(fc) 
and n{k) = n^ik) + ni{k) near R = Rc for PNNN-PIA 
and PNNN+PIA, respectively. They show that for a suf- 
ficiently small a, the degree of the congested node at 
R — Rc is generally close but not equal to fcmin. This 
is not surprising because there are numerous nodes with 
degree close to m. Local conditions such as the degrees of 
the neighbors of these small degree nodes can vary a lot. 
Combined with the break down of the scaling relation 
in Eq. pip , jamming may occur at a node whose degree 
is slightly greater than fcmin when R = Rc. In contrast, 
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rameters used are the same as those in Fig. 



Figs.[7|and[8]show that for a sufficiently large a, jamming 
almost always occurs in the highest degree node in the 
network. This is because for a generic BA network with 
a large but fixed N, there is a considerable difference be- 
tween the degree of the most connected and second most 
connected nodes. Thus, n; for the most connected node 
is almost surely greater than that for the slightly less 
connected ones. Most importantly, our simulations find 
that the transition between these two types of congested 
nodes aX R = Rc occurs at a rather well-defined critical 
ac for TV > 1000. And the value of etc depends on the 
value of N as well as on whether the PIA rule is adopted 
or not. 

After finish justifying the validity of the approxima- 
tions made in our mean field analysis, we now move on 
to compare our mean field calculations and numerical 
simulation results with the numerical findings of Yin et 
al. reported in Ref. [l3|. As the Rc against a curves in 
Figs. [5] and [TO] shown, the general trend of Rc we find 
in our numerical simulations agrees quite well with the 
predictions of our mean field theory for both PNNN-PIA 
and PNNN-I-PIA. In particular, we discover that for fixed 
A'^ and m, Rc is an increasing (decreasing) function of a 
for a < Uc {a > etc)- Besides, ac decreases and Rc in- 
creases with the introduction of PIA rule although the 



change is not significant for large N and small m. Re- 
call from Eq. (|30ap and the discussions in Sec. IIIIB II 
and IIII B] that in the large N limit, S should be roughly 
a constant over the parameter range of interest and Rc 
should roughly scales like 1/(1 — a) whenever a < ac- 
This is exactly what we find in Figs. [51 and [TOl More gen- 
erally, Eqs. pIa[) -([5Dc |) imply that RcN^'^'^^ should be TV 
independent, where 

{0 for a < -1, 

(a + l)/2 for -l<a< 1, (31) 
1 for a > 1. 

As shown in Figs, [ll] and [l2l RcN^^-"^^ is indeed TV inde- 
pendent for a < 1 (a > 1) provided that N > 10000 
(TV > 500000). Again, the discrepancy for a > 3 is 
probably caused by insufficient sampling and the finite 
precision of our pseudo random number generator. 

As for the critical preferential delivering exponent ac, 
we find that it decreases as m increases for a fixed 
TV. This can be explained as follows: Recall that 
the number of 4-cycles in a BA network scales like 
[mlog(TV)/2] /4 [1J|. So, by increasing m while fixing 
TV, the proportion of 4-cycles in the network increases. 
In other words, the assumption of neglecting the effect of 
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4-cycles in our mean-field analysis reported Sec. IIIII be- 
comes less valid. By going through the analysis in Sec. IIIII 
once more, it is not difficult to see that although the scal- 
ing relations in Eqs. pT|) and are robust against the 
presence of 4-cycles, the {k — 1) factor in Eq. should 
be replaced by {k — () for some C > 1- This change de- 
creases the value of ni(k) for a fixed N, therefore making 
the small degree node harder to jam. This is the reason 
why the presence of large number of 4-cycles reduces the 
value of Uc- 

In the case of to = 4, Figs. [D^a) and fTUf a) show that 
limjv^oo etc is very close to —1 for PNNNzbPIA. Com- 
bined with the validity of Eqs. (|30a[) and (j30bp as de- 
picted in Figs. [TlT a) and [121 we conclude that Rc ap- 
proaches to its maximum value at etc = —1 in the large 
N limit. In contrast, for simulation up to = 100000, ac 
does not seem to converge to —1 in the case of to = 5. As 
we have discussed in the last paragraph, we believe that 
this is due to the existence of large amount of 4-cycles. 
Since for m — 5, the number of 4-cycles is less than 
about A^/10 provided that A^ > 10'', we believe that ac 
should converge to —1 by using networks at least about 
100 times larger than our currently used ones. Unfortu- 
nately, such simulation is beyond the current computing 
capacity of our group. 

Now, let us compare our findings with that of Yin et 
aZ.'s in Ref. [T2|. Fig. [TUI show that the simulations per- 
formed on a A^ = 1000 network does not reveal the ther- 
modynamic behavior of the system due to serious finite 
size corrections. Actually, if they had extended their nu- 
merical simulations to a as small as about —8 (which un- 
fortunately requires much longer computational time and 
the use of a high precision pseudo random number gener- 
ator), they should have revealed the maximum point on 
the a — Rc curve, thereby discovering the critical ac- 

V. DISCUSSIONS 

To summarize, we have pointed out that the 
PNNN+PIA model is not a good model of network traf- 
fic due to the hidden communication cost involved. In 
addition, we have carefully performed a mean-field anal- 
ysis of the message packet dynamics for a network traffic 
model with PNNN routing strategy on BA network with 
or without PIA by Yin et al. in Ref. [l2|. The main fea- 
ture of our analysis is that we divide the message packets 
into two groups, namely, the DSPs and DLPs. To check 
the validity of our mean-field results, we introduce a new 
method to simulate the critical packet injection rate Rc 
that requires much less memory. This enable us to carry 



out an extensive numerical simulation to study the so- 
called a — Rc curve for larger network size A^ with the 
message packet injection rate R taking on real rather 
than integer values. 

For a fixed finite network size A^, we discover that 
the a — Rc curve is in fact increasing (decreasing) for 
a < ac {a > ac)- And we are able to explain this 
behavior by means of our mean- field analysis. In fact, 
both our mean-field calculations and our numerical sim- 
ulations show that the critical message generation rate Rc 
attains its maximum value at ac = — 1 for models both 
with and without PIA rule in the limit of large A^. In this 
respect, the role of PIA rule has little effect on the a — Rc 
curve even though the value of Rc is increased by intro- 
ducing the PIA rule. At the same time, Eq. (|30ap tells 
us that Rc is independent of A^ in the limits of A^ ^ oo 
and a —1^. This means that the PNNN mechanism 
is not efficient in handling large scale BA network traffic. 
In fact, this result agrees with those of Sreenivasan et 
al. [20] who showed that Rc < 0{\/N) for a BA network 
with any routing strategy. 

One may apply our analysis to consider the extension 
of PNNNzbPIA to the case in which more extended lo- 
cal information of the network such as the third nearest 
neighbors is used to forward a packet. It is not too dif- 
ficult to argue that ns(fc) ~ and n;(fc) ~ fc in the 
large A^ limit for this kind of models. Thus, it appears 
that straight-forward generalizations of the PNNN packet 
forwarding rule are also not efficient to handle large scale 
BA network traffic in the sense that the resultant maxi- 
mum possible value of Rc is independent of A^. One has 
to find other type of strategies in order to approach the 
upper bound of 0{^/N) for Re- 
in addition to the functional dependence of Rc on 
a, it is instructive to study the nature of the phase 
transition between the free-flow and jamming phases in 
PNNNzbPIA. Nonetheless, our mean field analysis and 
the trick used in our extensive numerical simulations are 
for free-fiow phase only. Further work has to be done to 
investigate this problem. 
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